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Abstract. Pulse walk-off in the process of sum frequency generation in a nonlinear 
•^(2) crystal is shown to be responsible for pulse jittering which is reminiscent to the 
Zitterbewegung (trembling motion) of a relativistic freely moving Dirac particle. An 
analytical expression for the pulse center of mass trajectory is derived in the no- 
pump-depletion limit, and numerical examples of Zitterbewegung are presented for sum 
frequency generation in periodically-poled lithium niobate. The proposed quantum- 
optical analogy indicates that frequency conversion in nonlinear optics could provide 
an experimentally accessible simulator of the Dirac equation. 
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1. Introduction 

Originally predicted by Schrodinger in the study of the Dirac equation pQ, 
Zitterbewegung (ZB) refers to the trembling motion of a freely-moving relativistic 
quantum particle that arises from the interference between the positive- and negative- 
energy parts of the spinor wave function [2j. For a free electron, the Dirac equation 
predicts the ZB to have an extremely small amplitude (of the order of the Compton 
wavelength ~ 10 -12 m) and an extremely high frequency (~ 10 21 Hz), making such 
an effect experimentally inaccessible. In addition, the physical relevance of ZB in 
relativistic quantum mechanics is a controversial issue because such an effect arises in 
the framework of the single-particle picture of the Dirac equation, but not in quantum 
field theory [31 0]. The notion of ZB and resulting formalism, however, are not peculiar 
to relativistic quantum dynamics, and phenomena analogous to ZB, which underly the 
same mathematical model of the Dirac equation, have so far predicted in a wide variety 
of quantum and even classical physical systems, including among others semiconductors 
and quantum wells [5j [6], [T] , trapped ions [8], graphene [91 UOj [Til [12] , cold atoms [T3irbf] . 
acoustic [15] and photonic [16j [T71 HE] systems. Simulations of relativistic quantum 
effects using experimentally-accessible physical set-ups, in which parameter tunability 
allows access to different physical regimes, have seen in recent years an increasing 
interest, culminating to the very recent first experimental observation of a quantum 
analogue of ZB using a single trapped ion set to behave as a free relativistic quantum 
particle ^L9\. In the optical context, the use of photonic systems to mimic quantum 
phenomena in the lab has seen a continuous and increasing interest (see, for instance, [20] 
and references therein); in particular, optical analogues of the relativistic ZB have been 
recently proposed to occur in photonic crystals [IB] , metamaterial slabs [17], and binary 
waveguide arrays [18] . In this work it is shown theoretically that a classical analogue 
of ZB can be observed in a much simpler and well-known set-up of nonlinear optics, 
namely in the process of sum frequency generation of light waves in a nonlinear x 
medium [21] in presence of pulse (or spatial) walk off. In the nonlinear optics context, 
optical three- wave interaction (TWI) in nonlinear x^ media in presence of temporal (or 
spatial) walk-off is a well-known process, which has been widely investigated especially 
in connection to pulse compression of ultrashort pulses and to TWI soliton theory (see, 
for instance, [221 1221 [2U [251 I2S1 [27] ) . Notably, the nonlinear TWI equations are solvable 
by inverse scattering methods [22] . However, the ZB phenomenon discussed in this work 
and the idea of exploiting nonlinear optics to mimic the Dirac equation have not been 
addressed in such previous studies. 

2. Basic model and quantum-optical analogy 

The starting point of our analysis is provided by a standard model of TWI in a nonlinear 
quadratic medium describing propagation of either optical pulses or optical beams at 
frequencies u\, 0J2 and 003 = u)\ + U2 in presence of either group velocity mismatch or 
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spatial walk-off. For the sake of definiteness, we will consider here the case of optical 
pulse interaction in presence of temporal walk-off (i.e. of group velocity mismatch), 
however the results hold also for spatial beam propagation in a quadratic medium 
with spatial walk-off provided that the temporal coordinate is replaced by a transverse 
spatial coordinate (see, for instance, [27]). Assuming that group velocity dispersion is 
negligible, in the plane- wave approximation and assuming perfect phase matching, pulse 
propagation in the x^ medium is described by the following set of nonlinear coupled 
equations (2Tj E21 E31 El EHl EE] : 

(> +~)A 1 = i P AiA a (1) 



dz v g i dt 

'' + —i) A s = ¥'AA (3) 



v dz v g 3 dt y 

where A\ = Ai(z,t) (/ = 1,2,3) is the amplitude of the electric field envelope at the 
carrier frequencies oji, normalized such that \Ai\ 2 is the photon flux at frequency ui, 
v g i is the group velocity in the medium at frequency U[, and p is the strength of the 
nonlinear interaction, which reads explicitly 



d e ff I 2fajJiU}2Ul3 



P=-^\ — ; (4 ) 

c V e c nin 2 n 3 

where n\ is the refractive index of the medium at frequency ui (/ = 1,2,3), Co is the 

(2) 

speed of light in vacuum, and <i e // = (l/2)x e ff is the effective nonlinear coefficient. To 
achieve perfect phase matching, a quasi-phase- matching (QPM) grating for the nonlinear 
susceptibility x can be employed; in this case one has (see, for instance, [28J) 

1- 



4// = -X (2) Wexp(zA^) (5) 

where Ak = (u^ris — U2U2 — uirii)/co is the phase mismatch of the three waves and the 
overbar denotes a spatial average over a few modulation periods of the QPM grating. 
Equations (1-3) admit of the following two invariants along the propagation distance z 

POO 

J 1>2 = / rft(|A li2 | 2 + |A 3 | 3 ) (6) 



which correspond to photon flux conservation (Manley-Rowe invariants) in the frequency 
conversion process. Solitary waves of Eqs.(l-3) in the fully nonlinear regime, including 
trapped bright-dark-bright solitary waves with locked velocity, have been investigated 
in Refs. [22l |23j EU EEJ ESI EZ]. To study the analogue of ZB in the frequency conversion 
process, we assume here that at the input plane z = the nonlinear crystal is excited 
by a strong and nearly continuous- wave pump field at frequency and by a weak and 
short signal pulse at frequency U2 and temporal profile g(t). Under such assumptions, 
the invariance of X 12 implies that the pump wave remains nearly undepleted along the 
propagation distance, and Eqs.(l-3) reduce to the following two linear coupled-field 
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equations describing sum-frequency generation in the undepleted regime 



fS 1 B\ , . , 



where we have set k = — pA\ = p^/ TjfajJi and I\ is the intensity of the pump field. 



Without loss of generality, k can be assumed to be real-valued and positive. Note that, 
in the absence of group velocity mismatch for the signal and second-harmonic waves at 
frequencies U2 and 0*3, i.e. for v 9 2 = v g3 , the solution to Eqs.(7) and (8) is analogous to 
the one for stationary fields [21], which shows a well-known oscillatory power transfer, 
along the propagation distance z, between the two fields with spatial period 7t/k; namely 
one has 



where v g = v g 2 = v g3 . Note that, in spite of the oscillatory power transfer in the 
frequency conversion process, the two pulses propagates with the common group velocity 
v g and do not show any trembling (jitter) motion. If the group velocity mismatch is not 
negligible (v g2 7^ v g3 ), the solution to Eqs.(7) and (8) is more involved, and is given by 
Eqs.(25) and (26) discussed in the next section. Here we anticipate that, in this regime, 
the oscillatory power transfer between the two fields is generally accompanied by an 
oscillatory motion of the pulse center of mass, which is reminiscent of ZB for the free 
relativistic Dirac electron. To highlight such an analogy in a formal way, it is worth 
introducing the coordinates of a moving frame 




(10) 



(9) 



f = z , rj = t- z/v g 
where the velocity v g is defined by the relation 




Vg 2 \Vg 2 Vg 3 J 

In the moving frame, Eqs.(7) and (8) take the form 




(12) 




(13) 



(14) 



where we have set 



2 \Vg 2 Vg 3 J 

Equations (13) and (14) are supplemented with the boundary conditions 
A 2 (o,v) = g(v) , A 3 (0, V ) = 0. 




(15) 



(16) 
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Figure 1. (color online) Sum frequency generation in a L = 1.5-cm-long PPLN crystal 
showing ZB of the signal pulse, (a) and (b): snapshots of the intensity distributions of 
signal [Fig. 1(a)] and sum-frequency [Fig. 1(b)] pulses versus propagation distance z. (c) 
Numerically-computed behavior of pulse center of mass (77)2 for the signal field versus 
propagation distance (solid curve), and corresponding behavior predicted by Eq.(28) 
(dotted curve, almost overlapped with the solid one). The inset depicts the behavior of 
the normalized photon fluence <\>i(z) / </>2(0) of signal held versus propagation distance, 
showing the oscillatory exchange of power between the signal and sum- frequency waves, 
(d) Numerically-computed behavior of (77), defined by Eq.(19), versus propagation 
distance (solid curve), and corresponding behavior predicted by Eq.(29) (dotted curve, 
almost overlapped with the solid one). The input pulse duration is t p = 5 ps. Other 
parameter values are given in the text. 



After introduction of the spinor wave field ip = (A 2 ,y4 3 ) T , Eqs.(13) and (14) can be 
finally cast into the Dirac form 

1— = -ia z b— + na x ip (17) 
a£ orj 

where a x and a z are the Pauli matrices. Note that, after the formal change 
8 -)■ c 

mc 2 , , 

(18) 

£ — > t , 7] —7- X 

Eq.(17) corresponds to the one-dimensional Dirac equation for a relativistic particle 
of mass m in absence of external fields, moving alon the x axis, written in the Weyl 
representation |2J. Therefore, the temporal evolution of the spinor wave function if} 
for the Dirac particle is mapped into the spatial evolution of the envelopes A 2 and A 3 
for signal and sum-frequency pulses, respectively, whereas the spatial coordinate of the 
Dirac particle is mapped into the retarded time 77 of the optical pulses. 
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Figure 2. (color online) Same as Fig.l, but for an input pulse duration t v = 1.5 ps. 




5 10 15 5 10 

propagation distance z (mm) 



5 10 

z (mm) 



5 10 

z (mm) 



Figure 3. (color online) Same as Fig.l, but for an input pulse duration r p = 0.5 ps. 
3. Zitterbewegung of optical pulses 

For the Dirac equation (17), ZB refers to the rapid oscillatory motion of the expectation 
value of the particle position 

n^(Mb| 2 + |A 3 | 2 ) 



(19) 



fZo dr )(\M 2 + IM 2 ) 

around its mean trajectory, which arises whenever negative- and positive-energy 
eigenstates of the Dirac equation are simultaneously excitated by the initial condition. 
Note that, indicating by (77)2(0 and (77)3(0 the (temporal) center of mass of the signal 
and sum-frequency pulses at the crystal plane £ = z, i.e. 

fZodnv\Awit,v)\ a 



one can write 



(vKt) 



02(0(^)2(0+03(0(^)3(0 



(20) 



(21) 



where 02,3(0 = / dr]\A 2 ,3(£,, rf)\ 2 are the photon fiuences of the signal and sum-frequency 
pulses at the plane z = £, respectively, and X 2 = 2 (O + 03(0 = 02(0) is the Manley- 
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Rowe invariant. In particular, if the initial pulse envelope g(t) is symmetric, i.e. 
g(—t) = g{t), as it will be shown below one has (77)3(0 = 0, and thus according 
to Eq.(21) the ZB of the Dirac particle for Eq.(17) can be simply retrieved from the 
temporal jitter (77)2(0 an d fractional energy 02(0/02(0) of the signal pulse. 
The solution to Eqs.(13) and (14) with 5^0 and with the boundary conditions (16) 
can be readily obtained in the spectral (Fourier) domain. Indicating by ^3(^,0;) = 
(l/27r) / dT/Aa^ (£, 77) exp(— irjco) the spectra of the signal and sum-frequency fields at 
the propagation plane £, one has 



A 2 (£,u>) = g{u>) 



cos(/30-iysin(/30 



(22) 



Mt,<")= -jg{w)tin(Pt) (23) 

where g(uu) = (l/27r) f dr]g(r]) exp(— ir/u) is the spectrum of the signal pulse incident 
onto the crystal at £ = 0, and 

p(u) = Vk 2 + u 2 5 2 . (24) 

In the temporal domain, the inverse Fourier transform of Eqs. (22) and (23) yields the 
following exact solution for the sum-frequency and signal pulses traveling along the 
crystal 



Si 
si 



MZ,v) = - l £l ^(0 + 77)Jo|M/e-^l (25) 



+ 



25 

9(V + 5£) + g(y - 5Q 
2 

I rSi d ( 





(26) 



where J is the zero-order Bessel function of first kind. To calculate (77), (77)2 and (77)3, 
let us assume, for the sake of simplicity, that the signal pulse envelope g(t) at the input 
crystal plane has a symmetric profile (with e.g. a Gaussian or a sech shape) satisfying the 
condition g(—rf) = g(rj). In this case, from Eq.(25) it follows that A 3 (£, —77) = A 3 (£, 77), 
and thus 

(^3(0=0, M(o=m<>7>2(o- (27) 

The explicit expression of (77)2(0; as obtained by substitution of Eq.(26) into Eq.(20), 
turns out to be rather cumbersome and not of easy physical interpretation. However, 
for a signal spectrum g(u) narrow at around 00 = with a spectral width Au much 
smaller than k/S, i.e. for a relatively long input pulse, simple approximate expressions 
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for (77)2(0 an d (t?) (0 can De obtained, which read explicitly 

5 3 Alu 2 



(vHO * £ 



cos 2 «) + 



<5 2 Aoj 2 



+ - Sln(2< L a (2* 
2k C os 2 (k0 + T 



fa>(0 2— + — sin ( 2 <) ( 29 ) 



k 2 Ik 

where Au is the spectral width of the input signal pulse, defined by 

9 f rfwa; 2 |9(a;)| 2 . , 

J du%(u;)| 2 

Equation (29) corresponds to the well-known approximate expression of ZB in relativis- 
ts quantum mechanics (see, for instance, [2J HH]), whereas Eq.(28) shows the signature 
of ZB in the oscillatory motion of the signal pulse center of mass as it propagates along 
the crystal. Note that such an oscillatory motion, that arises from the second term on 
the right hand side of Eq.(29), is superimposed to a slight drift term [the first term 
on the right hand side of Eq.(29)]. The oscillatory motion of (77) and (77)2 along the 
propagation coordinate £ of the crystal basically follows the oscillatory-like of optical 
transfer between signal and sum-frequency pulses. Note also that, according to Eq.(28) 
and because of the assumption 5Au/n <C 1, the pulse center of mass (77)2 (£ takes large 
values at the propagation distances £ = tt/2k, £ = tt/k, £ = 3tt/2k, ... where most of 
the signal field is converted into the sum-frequency field. 

As an example, let us consider the process of sum frequency generation in a nonlinear 
periodically-poled lithium-niobate (PPLNB) crystal (see, for instance, [29]) assuming 
Ai = 1550 nm, A2 = 810 nm, and A3 = 532 nm for the wavelengths (in vacuum) of 
pump, signal and sum- frequency waves, respectively. From Sellmeir equations [30], one 
can estimate at 25°C and for extraordinary waves ni = 2.1381, 722 = 2.1748, 773 = 2.2343, 
Vg2 /c = 0.4422, v g3 /c = 0.4069 and a QPM period of A = 2n/Ak ~ 7.39 /im, which 
is easily accessible with current poling technology. For first-order QPM with alternat- 
ing sign +/- of x with period A/2, the effective nonlinear coefficient is given by [29] 
d e ff = (2/7r)d, where d is the element of the nonlinear d-tensor of the crystal involved 
in the parametric interaction (d = d^ ~ 27 pm/V for extraordinary waves). As an 
input signal pulse, we assume a Gaussian profile g{t) oc exp(— t 2 /tq) with a full-width 
at half maximum (FWHM) pulse duration t p = (\/21og2)ro and spectral pulse width 
Au; = 1/tq. As an example, Figs. 1(a) and 1(b) show the evolution of the pulse intensity 
profiles |/l2(£,77)| 2 and |A 3 (£,7?)| 2 of signal and sum-frequency fields, respectively, in a 
L =1.5-cm-long PPLN crystal as obtained by direct numerical analysis of Eqs.(l-3), 
for a signal pulse duration t p = 5 ps of low peak intensity (1 W/cm 2 ) and an inten- 
sity of the continuous- wave pump field 1\ = ftwi\Ai\ 2 of 1 MW/cm 2 , corresponding 
to k — 0.4656 mm -1 and 8Au/k ~ 0.0827. The numerical results of Fig.l are with 
excellent accuracy reproduced by the analytical solutions (25) and (26), derived in the 
no-pump depletion limit. Figures 1(c) and 1(d) show the corresponding behavior, along 
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Figure 4. (color online) Numerically-computed behavior of the center of mass {rf) 2 
for the signal pulse at the output of the PPLN crystal versus the intensity I\ of the 
pump field. The crystal length is L = 1 cm, the input pulse duration of the signal 
field is t p — 5 ps. The dotted curve in the figure is the behavior of {rf}2 predicted by 
Eq.(28). 



the crystal coordinate £ = z, of the normalized photon fTuence 02(0/02(0) [inset of 
Fig. 1(c)], pulse center of mass (77)2(0 of signal field [solid curve in Fig. 1(c)], and ZB 
amplitude (77) (0 = [02 (0/^2(0)] (77)2(0 [solid curve in Fig.l(d)]. In Figs. 1(c) and (d), 
the behaviors of (77)2(0 an d ( r ?)(0> as predicted by Eqs.(28) and (29), are also shown 
(dotted curves, almost overlapped with the solid ones). Note that, according to the 
theoretical analysis, the center of mass of the signal pulse undergoes a clear oscillatory 
motion, superimposed to a slight drift [arising from the first term on the right hand side 
of Eq.(28)]. For spectrally broader signal pulses, ZB can not be accurately described 
by the simple Eqs.(28) and (29), however the oscillatory motion of the pulse center of 
mass, superimposed to a drift motion, is still observed in numerical simulations. This is 
shown, as an example, in Figs. 2 and 3, where pulse duration of the injected signal pulse 
has been reduced to t p = 1.5 ps in Fig.2 (corresponding to 5Au/k ~ 0.2758), and to 
t p = 0.5 ps in Fig.3 (corresponding to 5Au/k ~ 0.827). In an experiment, a measure- 
ment of the pulse center of mass at internal planes of the nonlinear crystal could be a 
nontrivial task, whereas autocorrelation measurements can readily reveal a jitter of the 
output pulse, at the exit of the crystal, with respect to a reference case. Owing to the 
dependence of the sinusoidal terms entering in Eq.(28) on the product k£ oc in an 

experiment the signature of ZB can be easier revealed by monitoring the center of mass 
of the signal pulse at the output plane £ = L of the crystal as a function of the pump 
intensity I\. This is shown, as an example, in Fig.4, where the numerically-computed 
behavior of (77)2 at the output crystal plane versus the intensity I\ of the pump field is 
depicted, together with the approximate prediction based on Eq.(28). 
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4. Conclusions 

In this work, a photonic analogue of the trembling motion (Zitterbewegung) of a free 
relativistic Dirac particle, based on frequency conversion of short optical pulses in a 
nonlinear quadratic medium, has been presented. The analogy, which stems from the 
mathematical similarity between the Dirac equation of a massive particle and the cou- 
pled equations describing sum frequency generation in the presence of pulse walk-off, in- 
dicates that well-known and experimentally accessible nonlinear optical processes could 
be exploited to simulate the Dirac equation in an optical setting. As compared to other 
classical and quantum analogues of Zitterbewegung recently proposed in the literature, 
based on trapped ions [19], graphene [9j ?,[TTJ[T2] or photonic crystals, superlattices 
or metamaterials [161 El EEE], our proposal may show a simpler experimental access 
and could stimulate further search for nonlinear optics analogues of relativistic quan- 
tum phenomena. For example, engineering of the QPM grating could be exploited to 
introduce in Eq.(17) a ^-dependence of k, i.e. to simulate the dynamics of a relativistic 
Dirac particle with a time- varying mass [31]. Likewise, if the temporal dependence of 
the pump pulse is included in the analysis and assuming v g \ = v g , an //-dependence of 
the mass k is introduced in Eq.(17), which enables to mimic a relativistic Dirac particle 
in a Lorentz scalar potential [2]. 

The author acknowledges financial support by the Italian MIUR (Grant No. PRIN-2008- 
YCAAK project "Analogie ottico-quantistiche in strutture fotoniche a guida d'onda"). 
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